Numerical simulations of two dimensional magnetic domain patterns 
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I show that a model for the interaction of magnetic domains that includes a short range ferro- 
magnetic and a long range dipolar anti-ferromagnetic interaction reproduces very well many char- 
acteristic features of two-dimensional magnetic domain patterns. In particular bubble and stripe 
phases are obtained, along with polygonal and labyrinthine morphologies. In addition, two puzzling 
phenomena, namely the so called 'memory effect' and the 'topological melting' observed experimen- 
tally are also qualitatively described. Very similar phenomenology is found in the case in which the 
model is changed to be represented by the Swift-Hohenberg equation driven by an external orienting 
field. 

PACS numbers: 



I. INTRODUCTION 

There is a surprisingly large number of systems that 
exhibit macroscopic textures arising from microscopic 
interactions.^ To be concrete, I will take as a case of 
study that of patterns in magnetic systems (magnetic 
garnets^' or ferrofluids"^), but many of the conclusions 
obtained can be directly applied to other systems, as 
for instance, the mixed state of type I superconductors 
of slab geometry)^ and Langmuir monolayers^ The phe- 
nomenology of these systems is qualitatively understood 
as appearing from the competition of two effects: a short 
range rigidity, and a long range (dipolar) interaction be- 
tween the local magnetization at different spatial posi- 
tions. Calculations suggest^ that the ground state of 
the system consists of (i) a state of uniform magneti- 
zation, (ii) a hexagonal lattice of bubbles in a back- 
ground with opposite magnetization, or (iii) a phase with 
alternating, parallel stripes of opposite magnetization. 
The parameter controlling which of these three is ac- 
tually the ground state is the external magnetic field. 
However, in experiments, upon variation of the external 
field, different (typically metastable) flux configuration 
develop that originate in instabilities of the bubbles or 
the stripes. Most noticeable, these metastable config- 
urations include labyrinthine phases of interpenetrating 
domains, and polygonal-like patterns.^ 

Model Hamiltonians that take into account the two 
relevant energy scales have been used to reproduce most 
of the elemental instabilities observed in experiments, in 
particular: the elongation and 'fingering' instability of 
bubbles,- and the undulation instability of stripes.^ How- 
ever, the much richer behavior of the full system, appear- 
ing from complex interaction effects in rather large spa- 
tial regions has not been studied in detail with this kind 
of models. In fact, it is not known if these simple mod- 
els contain all necessary ingredients to produce realistic 
magnetization patterns over large spatial scales. 

The main motivation of the present work is to present 
large scale simulations using a model Hamiltonian to see 
whether it can account for the full phenomenology and 



the variety of morphologies observed. I claim that the 
answer is positive. The simulations are able to reproduce, 
in particular, two phenomena that have been observed 
in these systems and have remained largely as puzzles, 
namely, the so called 'memory effect'— of some magnetic 
patterns, and the 'topological melting'^" of an ordered 
lattice of bubbles. 



II. DETAILS ON THE MODEL AND THE 
NUMERICAL TECHNIQUE 

The model I will use is not at all new (see [1], [11] 
and references therein). I will consider a scalar field (j){r) 
defined over the x-y plane. This variable will represent 
the magnetization in the system, that in experiments is 
typically constrained (because of structural properties) 
to point perpendicularly to the x-y plane. Then, exper- 
imentally, the magnetization (jj has a preference to take 
two different values, that without loose of generality I 
will assume to be ±1. It will be convenient for the simu- 
lations to consider (f) as a continuum variable and include 
in the Hamiltonian a local term Hi that favors the values 
6 — ±1. This term will be of the form 



Hi 



ao 



dr 



- ho / dr(/.(r) (1) 



and represents the simplest continuum field description 
of an Ising variable. Note that a term describing the 
effect of an external magnetic field ho has been already 
included. 

The other terms that will be included in the Hamilto- 
nian are the following. First, there is a rigidity term Hrig 
of the form 



Hr 



(3o / dr 



jV0(r) 



(2) 



This term (with positive /3o) discourages spatial varia- 
tions of (p, and can be called 'attractive', in the sense 
that two regions with a value of of plus (or minus) one. 
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in a background of the opposite sign, tend to merge into 
a single one to reduce the value of this term (in fact, in 
a description in terms of an Ising variable on a lattice, 
this term maps onto a ferromagnetic interaction between 
nearest neighbor sites). The fact that our fundamental 
variable (j) is continuous rather than discrete, and the 
existence of the gradient term, imply in particular the 
existence of a natural width (of the order of -y/po/o^) for 
the interface between domains with positive and negative 
magnetization. Choosing the parameters in such a way 
that this width is a few times the discretization distance 
in the simulation, allows to obtain a smooth interface be- 
tween domains, which turns out to be very weakly pinned 
by the underlying numerical mesh, and whose energy is 
almost independent of its spatial orientation. These two 
facts are crucial for a realistic simulation, and cannot be 
easily achieved using a Ising variable that takes only two 
values (see for instances the attempts in [12]).^'^ 

Secondly, there is a term H^^p that models the dipolar 
interactions, of the form 



Hd^p^lo / drdrV(r)0(r')G(r,r') 



(3) 



where G(r, r') ~ l/|r — r'|^ at long distances. At short 
distances however, the behavior has to be cut off to 
avoid divergences (in experiments, the cut off distance is 
given roughly by the thickness of the film). However, we 
can see that the way in which the cut off is done is not 
crucial for the results. In fact, we will take advantage of 
the fact that the two terms l(2Jl and (j3Jl can be compactly 
written in Fourier space as 



k 



|0(k)|2(/3ofc2+7oGk) 



(4) 



where Gk is the Fourier transform of G(r, 0). Thus, it is 
the combination (/3ofc^+7oGk) that will mostly determine 
the behavior of the system. Note that the short distance 
behavior of G in real space is masked in Fourier space at 
large k by the term, and then is irrelevant. On the 
other hand, the behavior at long distances transforms 
into a k dependence of the form 



G, 



The constants can be exactly evaluated to be 



/•OC 

uq ~ 271 rdrG{r) 
Jo 



fli = 27r 



(5) 

(6) 
(7) 



The finite value oq of Gk at fc = reflects the fact that 
the interaction in real space is integrable (in spite of be- 
ing sometimes called 'long range'). Also note that ai 
is independent of the short distance behavior of G(r). 
The main features of the interaction in Fourier space are 
the maximum with finite derivative at fc — > 0, and the 



minimum exists for any non-zero 70, indicating that the 
effect of the dipolar interactions on large distances can 
never be neglected. 

We have defined the energy function of the system, 
and now the dynamics has to be introduced. Since in 
magnetic systems the magnetization is a non-conserved 
order parameter, I will use the Allen-Cahn^** dynamical 
equations, namely 



dt 



(r) 



- -A (^ao(-0 + </-') -ho- (3oAcl) + -fo J dr>(r')G(|r - r'|( 

that represents an overdamped dynamics in which the 
system reduces its energy by a steepest descendant evo- 
lution. To efficiently implement these equations on the 
computer, and in order to avoid the direct evaluation 
of the integral in the last term of a pseudo-spectral 
methodic is used. I write the previous equation in Fourier 
space, namely 



9^ 
dt 



= -X [ao(-0 + 0')|k - ho6{k) + iPok^ + 7oGk)</)k] 

(9) 

In this way, the last term is now algebraic. The compli- 
cation has been translated to the first term, that involves 
the evaluation of the Fourier transform of (p'^. However, 
this can be done very efficiently by the use of standard 
fast- Fourier-transform techniques. 

In the simulations below, the function G is defined in 
real space to be G(r, r') = l/|r — r'p for any two points of 
the numerical mesh such that r 7^ r', whereas G(r, r) = 0. 
Then the cut off distance is the lattice discretization. 
The Fourier transform of this expression on the square 
lattice gives for the relevant terms of Gk the form of Eq. 
Q with Go — 9.05, ai = 2n. Once the value of oq is 
fixed (and since the value of ai is universal), there are 
four independent coefficients in l^. Two of them can be 
fixed by rescaling the spatial and temporal coordinates. 
In fact, if we define a new field 4>{r,t) = A~^(j){r/C,t/B) 
(and then (pk.{t) = A~^(j)ck{t/B j), and in case we choose 
A to be 



A = 



Qo7o(G - 1) 
ao 



(10) 



the new field satisfies equations of motion that in Fourier 
space can be written as (the tilde in the new field has 
been eliminated for simplicity): 



^ = a(0 - (/.3)|^ + M(k) - i(3k' + iGk)ct,^, (11) 



with 



minimum at a finite wave number fc„ 



7o//3o. This 



XuqA^ 
B 



(12) 
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(13) 
(14) 
(15) 



and where Gk is (up to the linear terms that are relevant 
for our analysis) the same function as before, namely 
Gk = ao — ai\k\ with the same gq and ai. This renor- 
malization can be used to fix two parameters in the new 
non-dimensional equations Hll|l . In the simulation pre- 
sented below I have fixed (3 = 2.0, 7 = 0.19 and took the 
spatial discretization to be the unit of length (this choice 
was convenient when implementing the equations on the 
numerical mesh, and have no other particular meaning). 
Therefore, we see that in addition to the external con- 
trol parameter h, a single internal control parameter a 
remains. This parameter regulates the possibility of the 
field (j) to take values others than the most convenient 
ones, namely 4> = ±1. We will see below the different 
morphologies that appear for different values of a. From 
now on I will always refer to the non-dimensional form 
(|11|1 of the equations of motion. 

Starting from an arbitrary initial condition, Eq. I|ll() 
describes an evolution in which the total energy of the 
system Hi+Hrig+Hdip is steadily reduced until it reaches 
a minimum, in which d4>k/dt is identically zero. We will 
see that typically the true minimum of the system is not 
reached, but instead one of many possible metastable 
states is obtained. The simulations presented below were 
done on a 512 x 512 mesh using periodic boundary con- 
ditions. The time-integration of the equations is done 
using a semi-implicit first order method, in which the fc^ 
term in Eq. (|ll|l is evaluated in the new time value. Con- 
cretely, I use an iteration scheme based on the following 
discretized form of (|ll|l 



St 



^'')\.+hS{k)~^Gu(t>i-f3k^4' 



(16) 

This treatment of the diffusive term is standard to im- 
prove the stability of the algorithm.i2. In all cases below 
the time interval used is St = 0.5. 



III. RESULTS 

The initial condition for the variable 4> is taken to be lo- 
cally random in the interval —1<4><1, and the system 
is evolved during an equilibration time tstart iu the pres- 
ence of a fixed applied external field hgtart ■ If h^tart is too 
large, the configuration obtained turns out to be a state 
of uniform magnetization. However, for lower hstart, a 
structure of bubbles of the minority phase (with magne- 
tization anti-parallel to the field) within a background of 
the opposite magnetization may be favored. 

After the time t start i the field is decreased as a func- 
tion of time with a finite rate dh/dt. This value is taken 




h= 0.008 



h=-0.0125 



FIG. 1: Evolution of the magnetization distribution upon 
reduction of the magnetic field h, for a = 1.6. Other param- 
eters are: tstart = 3000, tistart = 0.01, dh/dt = — 5 x 10~^ 
(see text). Here and in the following figures black (white) 
indicates regions with positive (negative) magnetization, all 
parameters are in the non-dimensional form corresponding to 
Eq. IIH . and system size is 512 x 512. 



to be as small as possible (within reasonable computing 
time) in order that the field change can be considered to 
be adiabatic (we will see that this cannot always be guar- 
anteed due to the existence in some cases of field driven 
instabilities). During the evolution, different morpholo- 
gies are observed for different values of a in Eq. Hll|l . 
which will be described now. 

Almost reversible interconversion of bubbles and 
stripes. For a = 1.6, the result obtained is shown in 
Fig. ^ Starting from the initial bubble phase, upon re- 
duction of the field /i, neighbor bubbles coalesce, forming 
a striped pattern. When the field becomes negative, the 
stripes destabilize, and separate in a chain of bubbles, 
which have opposite magnetization with respect to the 
original ones. The sequence of bubble and stripe pat- 
terns is found to be reversible upon cycling of the field. 



h=-0.0105 




h=-0.029 



h=-0.033 




h=-0.04 



h=-0.045 



h=0.03 



h=0.0075 




h=-0.024 



h=-0.0375 




h=-0.042 



h=-0.0465 



FIG. 2: Same as Fig. Qfor a = 1.8 {tstart = 1500, hs 
0.03, dh/dt = -3 X 10""). 



FIG. 3; Same as Fig. [Hfor a = 2.2 {t.tart = 4500, h, 
0.03, dh/dt = -1 X lO"**). 



There is however a noticeable hysteresis in the field value 
at which the bubble-stripe interconversion occurs. This 
is just the consequence of the transition between bubbles 
and stripes being first order. ^ 

From bubbles to rather isolated and wandering stripes. 
For a slightly larger value of a, namely a = 1.8 (Fig. 
O, the bubbles may become unstable and elongate in- 
dividually, without merging with their neighbors at the 
beginning. When they finally merge (for h ^ —0.025), re- 
gions with positive magnetization generate wavy stripes 
of well defined thickness. Contrary to the previous case, 
these regions do not separate into 'beads' when the field 
is made more negative, but eventually retract back to 
a single spot of positive magnetization that eventually 
disappears. 

Collapse of the bubbles to a polygonal pattern. For 
larger values of a, the bubbles are seen to remain (meta-) 
stable down to a field where they start to merge with their 
neighbors, but now in a sort of two dimensional way, as 
seen in Fig. Ofor a = 2.2. This has to be compared with 
the previous case where the initial collapse of bubbles was 



mainly one-dimensional, generating stripes (see the cases 
h = -0.029 and h = -0.033 in Fig. EJ. In the present 
case, the collapse of neighbor bubbles seems to occur as 
a cascade process, where some initial coalescences trigger 
the full transition of the lattice. In fact, in Fig. 0]the 
field was kept constant at the value h = —0.0465 (corre- 
sponding to the last panel in Fig. and the evolution 
was followed as a function of time. A coarsening process 
is occurring here. Actually, the last pattern in Fig. ^is 
not totally relaxed yet. Incidentally note in the last panel 
of Fig. ^the existence of small pentagonal bubbles, high- 
lighted by the arrows. This structure has been observed 
experimentally to be ubiquitous, and very stable. 

This case suggests the following interesting result: if 
a perfect original pattern of bubbles is constructed by 
hand, it can remain stable for values of the field at which 
the disordered bubble system would have already col- 
lapsed. Now, if in this ordered, metastable structure, a 
defect is introduced, it can completely disorder the lat- 
tice. In fact, we see in Fig. |Slhow the presence of the 
defect produces a sequence of instabilities that destroy 



t=2250 t=6750 




t=27000 t=l 10000 

FIG. 4: The final configuration if Fig. |21 evolved at con- 
stant field h — —0.0465 as a function of time, as indicated 
{t — corresponds to the last panel in Fig. |5Jl. Arrows in the 
last panel highlight some small pentagonal bubbles, a struc- 
ture that appears ubiquitous both in experiments an in the 
simulations. 



many of the walls between neighbor bubbles, generating 
a rather well defined disordering front that leaves behind 
a disordered structure with much lower magnetization. 
This effect has been experimentally observed and called 
topological melting^ of the bubble lattice. It has been ob- 
served to occur (although in a less dramatic form) also for 
systems in which the long range interaction is of Coulomb 
typeii. 

Labyrinthine patterns and the memory effect. If from 
the last panels in Figs. |21or01the field is slowly switched 
off, interesting results are obtained. In the case in which 
we start from the configuration of the last panel in Fig. 12 
which contains three (meta-) stable spots of positive mag- 
netization, they remain stable (increasing only slightly in 
size) up to /i ~ —0.02. At this field an instability occurs, 
the bubbles becoming unstable. If we maintain the field 
fixed at a value slightly lower (in absolute value) than the 
instability value, we obtain the results presented in Fig. 
El which shows snapshots as a function of time, for a fixed 
value of the field h = —0.018. The bubbles elongate and 
successively branch, forming labyrinthine patterns that 
invade the whole sample. On the other hand, if the field 
that we apply is much beyond the instability value (Fig. 
I?)) the evolution is more rapid, and with a larger degree 
of branching of the magnetic domains. Note the differ- 
ence in the degree of branching in the final patterns of 
Figs. El and The grater tendency to branching when 



t=0 t=1500 




t-5500 t=9000 




t=13000 t=25000 

FIG. 5: Topological melting of an order array of bubbles for 
a = 2.2, upon the ad hoc inclusion of a defect in the middle of 
the sample. The evolution occurs at a fixed value of the field 
h = —0.05, as a function of the simulation time, as indicated. 



the applied field is more and more beyond the instabil- 
ity value is well known experimentally and theoretically^. 
This kind of instability is also similar to that observed is 
some reaction-diffusion systemsASi. 

If we reduce the absolute value of the field from a con- 
figuration in which stripes are already present, we observe 
an undulation transition^ at a field with larger absolute 
value than before. But contrary to what happened in 
Figs El and [3 if the field is changed slowly the system 
evolves smoothly (no instability appears), and stripes 
do not branch. Positive magnetization regions invade 
the system through wandering of the stripes, but new 
branches do not appear or are very rare. 

In particular, in the case in which we reduce the field 
starting from the last configuration in Fig. 01 in which 
stripes are abundant as walls between polygons, the un- 
dulation occurs mildly, with almost no breaking or recon- 
nection of the cell walls, and then the final labyrinthine 
pattern at /i = is topologically equivalent to the orig- 
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FIG. 6: Time evolution of the pattern shown in the first panel 
upon the application of a constant field {h — —0.018) slightly 
beyond the critical field at which those bubbles destabilize. 
Times of the snapshots are indicated. 



FIG. 7: Same as Fig. |S]for h = 0, i.e., here the system is 
brought deeply inside the instability region. Note the much 
larger amount of stripe branching in the final (stable) pattern, 
and the shorter time scale as compared with the previous 
figure. 



inal one. This is shown in Fig. |S1 A nice consequence 
of this, is that when the field is switched on again (last 
panels in Fig. IHJ, the original pattern is almost recov- 
ered. This effect, called the memory effect has been 
observed experimentally, and the typical evolution of the 
patterns over many cycles of the field has been analyzed. 
We are seen that this effect is contained in the simple 
model Hamiltonian we are using. 

IV. DISCUSSION AND CONCLUSIONS 

Summarizing, in the previous Section I have shown 
how the model equations can be efficiently simu- 
lated in systems of reasonably large size. In this way, we 
have seen emerging most of the phenomenology of two di- 
mensional magnetic patterns and other similar systems. 
The success of the present numerical simulations are due 
to a combination of reasons, mainly: the use of a contin- 
uum variable instead of a discrete one to obtain smooth 
domain walls between regions with opposite values of 0, 
and the use of pseudo-spectral techniques to evaluate ef- 
ficiently the 'long-range' dipolar force. These facts com- 
bine to allow a realistic simulation of domain patterns 
that show many of the features observed in experimental 
realizations. In particular, the memory effect and the 
topological melting^^ of the system are very well repro- 
duced. 

I want to emphasize that in all cases I have studied, the 



evaluation of the total energy of the system is compatible 
with the fact that the only patterns truly corresponding 
to the ground state of the system are: (i) a pattern with 
uniform magnetization if the field is strong enough, (ii) 
a regular bubble phase for intermediate fields, and (iii) 
a regular stripe phase for low (including zero) field. Al- 
though this is not a demonstration that they are the only 
possible ground states, it points in this direction, and it is 
in agreement with the results of theoretical studies.® The 
other patterns observed (labyrinthine, polygonal, etc.) 
are seen to be metastable, and they are originated in the 
particular cycling of the field (and in the initial condi- 
tions) to which the sample is subjected. A recent exper- 
imental study^^ has shown in fact how the labyrinthine 
patterns converge to parallel stripes upon relaxation. 

Very different morphologies have been observed when 
the parameter a in Eq. l(TT|l is changed. Figs. 01 S 
and |S1 (corresponding to the largest values of a) compare 
very well with the patterns observed in magnetic garnets 
and ferrofluids (see [1], [10], and [11]). The results for 
lower a (in particular. Figs. Q and |2Jl are more akin 
to Langmuir monolayers-"? and flux structures in type I 
superconductors.^ This suggests that in real systems the 
possibility of the order parameter to take values different 
than the two preferred ones can influence noticeable the 
physical properties. 

I want to mention here that the present model can be 
also efficiently used to study the effect of quenched disor- 
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h=-.0465 h=-0.0285 




h=-0.0315 ]i=-0.0465 

FIG. 8: Reducing the field from the final configuration in 
Fig. 0|down to /i = and back to its original value {\dh/dt\ — 
1 X 10~®). Note the 'memory' of the pattern as comparing 
first and last panels. 

der in the system, and the effect of thermal fluctuations. 
Preliminar results indicate that the model generates hys- 
teresis curves and magnetization patterns that, as a func- 
tion of the amount of disorder, compare very well with 
experimental ones^S. These results will be published sep- 
arately. 

We have seen that in the present model the dynamics is 
controlled by an interaction function in k space that has 
a maximum with finite derivative at fc — > and a mini- 
mum at a finite fc,„i„ value. It is worth comparing this 



case with respect to other possibilities. One is the case 
in which the field cj) is considered to be charged, instead 
of carrying a dipole. Two cases can be considered. One 
is that of true three dimensional charges (G'(r) ~ r~^) 
and the other is the case of two dimensional charges 
(G(r) ^ — ln(r)). In both cases, the interaction in k 
space gets a divergence at low k. This model has been 
studied in detail in [17] (see the references in there for 
realizations of this case). There, instabilities of a sin- 
gle bubble have been found which are similar to those I 
find in the dipolar system. It remains to be seen if the 
other effects described here are also present in Coulombic 
systems. 

Another case to compare with is that of interactions 
decaying in real space more rapidly than r^^. In this 
case, a k space interaction with a quadratic maximum 
at fc = is obtained. If this maximum dominates over 
the quadratic minimum coming from the A(j) term in 
Eq. (|SJ), then the effective interaction has a quadratic 
maximum at the origin and a minimum at some finite 
kmin- This case corresponds qualitatively to the interac- 
tion considered in the Swift-Hohenberg equation. For 
this interaction, and controlling the same parameter a 
as I did here, I have obtained basically all the effects 
and morphologies described in the previous section. On 
one hand this tells that the singularity at fc = of the 
dipolar interaction is not crucial in obtaining these ef- 
fects, a quadratic maximum suffices. On the other hand 
it is a bit surprising that in the wide literature related 
to the Swift-Hohenberg equation these effects have not 
been described previously. This might be due to the fact 
that the Swift-Hohenberg equation is usually considered 
in the absence of a 'magnetic-field-like' term that favors 
one of the two orientations, and this term is crucial to 
obtain the metastable patterns. It is then likely that the 
much studied relaxation to equilibrium properties of the 
patterns seen in the Swift-Hohenberg equation and the 
coarsening properties of the magnetic patterns (studied 
for instance in [9]) can be put under the same framework. 
I hope the present work encourages some studies in this 
direction. 
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